%% stagedOptimization_v1_3
%  Version 1.3
%  Author: Adeyinka Lesi
%  Date: 3/17/20
%  Project: Parameter optimization
%  Description: run optimization at different jump intervals and multiple
%  initial points to try to produce a global minimum
%% Version History
%  1.1: adjusting to bremermann_parallel_v2_0, which has input and output
%  history
%  1.2: more rounds - exploring more areas of the map
%  1.3: include scaling factor for levs
%  4/21/22: changed ds prefactor from 0.25 to 0.10

function [x1,fval,flag,output] = stagedOptimization_v1_3(opt_func,func,x0,low,high,scale_levs,opts)

start = clock();
% set defaults
if(~exist('opts','var'))
    opts = struct('del',0.05,'maxIter',1000,'TimeLimit',24*3600,...
        'Tolerance',1e-8,'ToleranceLength',10);
end
if(~exist('scale_levs','var'))
    levs = 5;
    scale = 1.12;
else
    levs = scale_levs(2);
    scale = scale_levs(1);
end
output = struct('iterations',0,'message','',...
    'totaltime',0,'bounds',[low,high],...
    'history',[],'inputs',[],'outputs',[]);

N = opts.maxIter;
ds = 0.10*scale.^(1-(1:levs));
n = floor(N/levs);
rem = N/levs-n;
if(rem ~= 0)
    nrem = abs((rem)*(1:levs)-round((rem)*(1:levs)))<1e-10;
    dN = ones(1,levs)*n+nrem;
else
    dN = ones(1,levs)*n;
end
% range = high-low;

x1 = x0;
fval = inf;
npas = 0;
history = zeros(1,N+levs);
outputs = zeros(1,8*N+levs); % assuming bremermann with 5 points
inputs = zeros(size(outputs,2),length(x0));
neval = 0;
for lv = 1:levs
    optlv = opts;
    optlv.del = ds(lv);
    optlv.maxIter = dN(lv);
    optlv.TimeLimit = opts.TimeLimit/levs;
    [xlv,flv,flag,outlv] = opt_func(func,x1,low,high,optlv);
    if(flv <= fval)
        x1 = xlv;
        fval = flv;
    end
    history(npas+lv:npas+lv+outlv.iterations) = outlv.history;
    npas = npas+outlv.iterations;
    if(isfield(outlv,'outputs'))
        nevallv = length(outlv.outputs);
        outputs(neval+(1:nevallv)) = outlv.outputs;
        inputs(neval+(1:nevallv),:) = outlv.inputs;
        neval = neval+nevallv;
    end
end

output.iterations = npas;
output.message = [num2str(levs) ' levels: ' outlv.message];
output.history = history(1:npas+levs);
output.totaltime = etime(clock(),start);
output.inputs = inputs(1:neval,:);
output.outputs = outputs(1:neval);